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ABSTRACT 


A  method  is  presented  for  obtaining  numerical  solutions  of  the 
Eulerian  hydrodynamic  equations  for  supersonic  compressible  flows.  The 
procedure  is  an  adaptation  of  the  PIC  method  to  the  case  of  a  continuous 
fluid.  Some  discussion  is  given  concerning  the  limitations  of  the  pro¬ 
cedure;  the  conservation  of  mass,  energy,  and  momentum;  and  the  treatment 
of  rigid  obstacles  and  two-fluid  interfaces.  Examples  of  calculations 
for  flows  past  a  number  of  cylindrically  symmetric  objects  are  also  pre¬ 
sented. 
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I.  INTRODUCTION 


This  report  has  been  written  with  a  two-fold  purpose.  The  first  is 
to  record  a  set  of  procedures  developed  several  years  ago  for  the  numeri¬ 
cal  solution  of  the  Eulerian  hydrodynamic  equations  for  compressible  isen- 
tropic  flow.  The  second  is  to  present  the  results  of  some  recent  computa¬ 
tions  on  the  supersonic  flow  of  air  past  a  variety  of  cylindrically 
symmetric  objects. 

The  procedure  for  solving  the  hydrodynamic  equations  to  be  described 
here  is  closely  related  to  the  semi-Eulerian  "Particle-in-Cell"  (PIC) 
method1  and  is  in  fact  merely  an  extension  of  that  approach  to  a  contin¬ 
uous  fluid.  Like  the  PIC  method,  it  is  applicable  to  a  wide  variety  of 
multi-fluid  hydrodynamic  problems  in  one  or  more  dimensions  and  is  par¬ 
ticularly  suited  to  situations  involving  large  distortions  of  the  fluid. 
Limitations  of  the  procedure  are  imposed  primarily  by  bounded  instabili¬ 
ties  which  arise  when  large  portions  of  the  fluid  under  consideration 
move  at  velocities  below  the  local  sound  speeds.  Otherwise,  calculational 
accuracy  is  governed  mainly  by  the  fineness  of  the  Eulerian  difference 
mesh  one  must  use,  which  is  determined  by  the  storage  capacity  and  speed 
of  available  digital  computers. 
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Sections  II  and  III  of  this  report  are  devoted  to  a  discussion  of  the 
numerical  techniques  involved  in  the  solution  of  the  hydrodynamic  equa¬ 
tions.  In  Section  IV,  theoretical  results  are  given  for  steady- state  flow 
around  right  circular  cylinders  with  conical  noses  of  various  apex  angles 
and  at  various  Mach  numbers  between  one  and  five.  Comparisons  of  the  com¬ 
puted  positions  of  detached  shocks  and  of  the  pressure  distribution  about 

2 

the  obstacle  nose  are  made  with  the  experimental  results  of  Marschner  and 
3 

Oguchi. 


II.  CAlCULATIOmL  PROCEDURE 

A.  Description  of  the  Method 

The  method  for  solving  the  hydrodynamic  equations  which  is  described 
below  finds  its  primary  applications  in  two-dimensional  problems.  We  will 
therefore  forego  a  discussion  of  the  procedure  in  one  dimension  and  begin 
immediately  with  a  discussion  of  the  two-dimensional  differencing  scheme. 
For  simplicity,  only  a  single  fluid  in  Cartesian  coordinates  will  be  con¬ 
sidered  at  this  time. 

In  order  to  obtain  a  system  of  difference  equations  to  describe  the 
motion  of  the  fluid,  we  assume  that  the  region  of  fluid  under  considera¬ 
tion  is  covered  by  a  fixed  (Eulerian)  grid  of  rectangular  cells  of  sides 
A x  and  Ay  (Fig.  1 ). 


* 
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Figure  1 •  Example  of  a  tiro  dimensional  Eulerian  mesh 

Each  cell  is  labeled  by  the  pair  of  integers  (i,j)  which  denote  the  loca¬ 
tion  of  the  cell  center  with  respect  to  some  origin.  Half  integer  indices 
are  used  to  denote  cell  boundaries.  For  example,  the  indices  ( ±+it 0 )  re_ 
fer  to  the  right-hand  boundary  of  cell  (i,j),  and  (i, j-J)  refers  to  the 
lower  boundary  of  this  cell.  A  calculation  is  begun  by  assigning  values 
for  the  density  p,  velocity  v*,  and  specific  energy  E  of  the  fluid  at  the 
center  of  each  cell  at  the  initial  time  t  =  0.  From  these  initial  condi¬ 
tions  the  configuration  of  the  fluid  over  the  mesh  is  obtained  at  a 
slightly  later  time  t  =  At  by  means  of  the  hydrodynamic  equations  of  motion. 
In  general,  one  proceeds  cyclically  to  calculate  the  state  of  the  fluid  at 
time  t  =  nAt  from  the  state  at  time  t  =  (n-1  )At  where  n  is  an  integer 
which  will  henceforth  be  used  without  the  At  to  denote  the  time. 

In  integral  form,  the  hydrodynamic  equations  for  compressible  isen- 


tropic  flow  are 


a 

/  pv  dT  =  - 

Jv(t) 

f  p 

s(t) 

(2) 

* 

a 

3t 

/  E  dT  =  - 

'v(t) 

f  p  V  •  ds 

7S(t) 

(3) 

p  = 

p(p,l) 

j 

w 

Here,  we  are  considering  integration  over  a  specific  element  of  fluid  oc¬ 
cupying  volume  v(t)  at  time  t  and  bounded  by  surface  S(t).  It  is  assumed 
that  no  external  forces  or  internal  dissipative  mechanisms  are  present 
anci  that  the  motion  of  an  element  of  fluid  is  produced  solely  by  the  pres¬ 
sure  p  exerted  on  this  volume  by  neighboring  fluid  elements.  The  first 
three  equations  are  statements  of  the  conservation  of  mass,  momentum,  and 
energy  of  the  fluid;  and  the  last,  the  equation  of  state,  relates  the 
pressure  of  the  fluid  to  the  density  and  the  specific  internal  energy, 

I  =  E  -  i  v2. 

In  differencing  these  equations,  we  will  consider  the  volume  V(t) 
to  be  that  of  some  arbitrary  cell  (i,j)  and  the  surface  S(t)  to  represent 
the  four  sides  of  this  cell. 

As  in  the  PIC  method,  the  calculation  of  the  time  (n+1)  configura¬ 
tion  of  the  fluid  from  that  at  time  n  may  be  thought  of  as  proceeding  in 
three  steps.  In  the  first,  the  fluid  element  in  each  cell  is  assumed 
fixed  and  no  flow  across  the  cell  boundaries  is  allowed.  By  applying 

Eqs.  (2),  (5 ),  and  (4),  intermediate  values  of  the  velocity  and 
specific  energy  in  each  cell  are  obtained.  The  second  phase  consists 
of  the  calculation  of  the  mass  flows  across  the  cell  boundaries. 
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Mass  crossing  a  boundary  is  assumed  to  carry  with  it  the  intermediate 
velocity  and  specific  energy  of  the  cell  from  which  it  came.  In  the 
third  phase,  the  new  cell  densities  are  computed;  and  by  application  of 
conservation  of  energy  and  momentum,  the  new  velocities  and  specific 
energies  are  subsequently  obtained. 

In  the  following,  all  cells  under  consideration  will  be  assumed 
interior  to  the  fluid.  The  effect  of  boundaries  will  be  left  for  later 
discussion. 


B.  The  Difference  Equations 


1 .  Phase  I 

It  is  assumed  that  the  following  quantities  are  known  at  the  center 
of  each  cell  (i,j)  at  time  n,  i.e.,  t  =  nAt,  and  that  for  purposes  of 
differencing  they  remain  constant  throughout  the  cell: 


density: 

U) 

velocity: 

u(n> 

M 

(x  component) 

v<a>. 

(y  component) 

specific  energy: 

(a) 

By  applying  the  equation  of  state  Eq.  (4),  the  pressure  at  the 
center  of  each  cell  is  obtained. 


=  P 


2 

2 


fn)£ 

Lui,j 


21 


+  v 


(n) 
i,j  J 
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Integrating  Eqs.  (2)  and  (3)  over  the  volume  of  cell  (i,j)  at  time  n  re¬ 


sults  in  the  difference  equations. 


(n)  /duN(n)  A  A  (n) 

pi,J  (3t>'  AxAy  '  -  Lpi4,j 

J 


P<nl 

1-2, JJ 


Ay 


(n)  /5vN(n)  .  . 

4J  (stK  ’ AxAy  = 


(n)  _  (n) 

Pi,j+|  Pi,j-i 


Ax 


(n)  /^E\(n) 

*’J  *1,3 


AxA  y  =  - 


Jn) 
P-t  j_I 


(n) 
u.  i 


>)  >> 

”  P-J  1  A  U„.  1 


l  i+i, j  ±+2,j  *-2>J 


Ay 


(n)  (n)  (n)  (n)  . 

v)  ',1  v'  ;  i  -  v  '  1  v)  '  1 

-  2  1>J+2  i#J“2  ■'■jJ  2J 


Ax 


The  assumption  of  no  mass  flow  has  been  used  to  remove  the  density  from 
the  time  derivative  on  the  left-hand  sides  of  these  equations.  The  pres¬ 
sures  and  velocities  on  the  right  refer  to  the  cell  boundaries  and  are 
taken  as  simple  averages  of  these  quantities  in  the  cells  separated  by  a 
given  boundary.  For  instance. 


1 


Pi+i,j  2 


(  n)  ( n) 

,d'  '  +  pv  • 

LPi,j  Pi+1,j. 


Defining  the  time  derivatives  on  the  left  sides  of  these  equations  by 


/  .  ~(n)  (n) 

-v  (n)  u.  .  -  u:  ' 

(|)  . 


At 


.  (n)  v(n)  -  v(n) 

&)  =  A/.J . i,.J. 

*  i,J 


At 


0 
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we  arrive  at  the  following  equations  for  the  intermediate  velocity  and 


specific  energy,  which  are  distinguished  by  a  tilde  symbol: 


v. 


E. 


-  u(n)  -  ' 

(n)  _  (n) 

‘i+Lj  Pi-^0  At 

(5) 

Ax  Tn) 

-  >) . ! 

PH  i  -  pH  i  A+ 

i» J+?  At_ 

(6) 

Ay  IS) 

M 

rPH,  .  uH  .  -  pH  •  pH  • 

1+PfJ  l+pfj  l"?fj 

-  - 

"  Ei,j 

Ax 

+ 

PH.  1  vH  !  -  PH  !  vH.  n 

tfJ+2  lfJ+2  hJlzz _ ~klA~-Z 

Ay  J 

At 

(7) 

2.  Phase  II 

From  the  right-hand  side  of  Eq.  (1)  it  is  seen  that  the  flow  across 
a  boundary  from  one  cell  into  another  is  proportional  to  the  fluid  den¬ 
sity  and  normal  component  of  velocity  at  the  interface.  Equation  (l) 
differences  into  the  following  equation  expressing  the  conservation  of 


mass: 


AxAy  =  p<“>  AxAy  -  ♦  A^J  -  A*^  +  AM«4 


(8) 


where 


(9) 


A„(n)  (n)  ~(n)  .  A. 

AM  I  .  =  pVi  .  u.  {  AyAt 


AMi”U 


p(n}  1  1  Ax  At 

i,0+j>  i,J+S 


etc.  are  the  mass  flows  across  the  cell  boundaries  indicated  by  the  sub¬ 
scripts. 

It  has  been  found  that  it  is  now  no  longer  adequate  to  use  simple 
averages  for  the  densities  and  velocities  at  the  cell  boundaries  as  was 
done  in  Phase  I.  Doing  so  tends  to  lead  to  minor  but  noticeable  insta¬ 
bilities  in  the  solution  of  the  difference  equations.  Instead,  a  weight¬ 
ing  must  be  given  these  quantities  in  preference  of  the  cell  out  of  which 
the  fluid  flows.  (This  necessity  arises  "naturally"  when  one  remembers 
that  due  to  the  finite  time  step.  At,  the  center  of  the  fluid  element 
leaving  a  cell  lies  somewhat  within  the  cell  and  not  at  the  boundary.) 

In  the  mass  flow  calculation,  the  boundary  quantities  can  be  obtained 
from  the  first  two  terms  of  a  Taylor  expansion  about  the  center  of  the 
emptying  cell.  Thus  the  equations  for  the  boundary  quantities  depend 

on  the  direction  of  flow  at  the  cell  interface.  For  the  right-hand  boun- 

/  \ 

dary  of  cell  (i,j),  for  instance,  the  mass  flow  AM.  {  .  would  be  calcu- 

1+2,0 

lated  as  follows: 


If 


.  ~(n) 

+  u 

i+1,j 


>  0  and 


~( n)  ~( n) 
uj  .  i  .  =  uj  ' 


1+2,  0 


i,0 


+  W.  . 

i,0 


Ax 

2 


~( n) 

=  u)  + 

1,0 


r(n) 


u:  .  . 

X+T 1 0 


u 


(n) 


>  0 


then 
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AM^}  ,  , 

1+2,  j  L  i >3 


,^n)  _  ^n) 


(n)  _  Jn) 


~(n)  +  VlfJ  “  Vl,jl  (n)  +  Pi+1  ,g  ;  Pl-M 


J  LHi,j 


AyAt 


If  ^n),  +  u^.n,  .  <  0  and 


+Vi,j 


(n)  <v  (a) 

-  s<“l  ,  -  (^)  lt!H  M<° 


i+|,0  i+1,0 


i+1,J 


then 


(n) 


*«i:L  -  -  V2'\'  ^ 


r(n).  r  (n)  _  (n). 

(n)  .  pi+2,j  Pi,j 

Lpi+1 , j  4 


AyAt 


Otherwise  AM^I  .  =  0. 

1+2>0 


The  flows  across  the  other  three  sides  of  cell  (i,j)  are  obtained  in 
a  similar  manner.  The  requirement  that  the  sign  of  the  simple  average 


boundary  velocity 


,u^n^  +u|Dl  in  the  case  above 

L  i,j  i+1, j 


be  used  in  determin¬ 


ing  the  form  for  the  mass  flow  term  ensures  that  a  calculation  will  be  in¬ 
dependent  of  the  direction  one  chooses  in  progressing  from  cell  to  cell 
through  the  mesh. 

When  more  than  one  fluid  is  present  it  is  convenient  to  drop  the  de¬ 
rivative  part  of  the  density  term  in  the  above  mass  flow  equation,  leaving 
only  the  first  term  in  the  Taylor  expansion.  This  does  not  noticeably 
affect  the  stability  of  the  difference  equations. 


3.  Phase  III 

Mass  flowing  from  one  cell  to  another  is  assumed  to  carry  with  it  the 
tilde  velocity  and  specific  energy  of  the  cell  from  which  it  came.  Thus, 
after  the  mass  flow,  the  composition  of  a  given  cell  may  be  thought  of  as 
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having  a  fraction  with  mass  at  velocity  v1  and  specific  energy  E^ ,  a 
portion  U  with  velocity  v^  and  specific  energy  E^  etc.  A  single  ve¬ 
locity  and  specific  energy  for  the  cell  can  he  obtained  by  applying  the 
laws  of  conservation  of  moment-urn  and  energy: 


To  illustrate  this  by  a  specific  example,  consider  the  changes  in 
cell  (i,j)  due  to  the  transport  of  fluid*  We  assume  that  the  fluid  flow 
in  the  neighborhood  of  this  cell  is  upward  and  to  the  right*  Therefore, 
fluid  flows  into  (i,j)  across  its  lower  and  left  boundaries  and  flows  out 
across  its  upper  and  right-hand  sides*  The  composition  of  the  cell  after 
transport  is  as  follows: 


mass: 


M1  =  AM^l  , 
1  1-2, J 


Mp  s  AM  \n\  i 


M,  = 


3  L  i,  j 


pfn^.  A  xAy  -  AM^l 


1+2,  J 


Al|n^.  i 


y 


* 


-1 6- 


carrying  velocity: 


?r(n)  ~(n) 

L  i-1  ,y  Vl,J. 


u: 


(n)  ~(n) 

i,  j-1 ’  i,j-lJ 


u 


(n)  ~(n) 


L  1,0 


i,0J 


carrying  specific  energy: 


■r(n) 


E 


(n) 

i^O-1 


E^n). 

i>0 


The  time  (n+1)  parameters  for  cell  (i,j)  will  thus  be 


(n+1) 

Pi,d 


(n+1) 
u.  . 


(n+1) 
v.  . 


M.  +  M-  +  M_ 

_  1  2  2 

Ax  Ay 

M  u(n]  +  Mu?n?  +  M  u^n? 

=  1  i-1 , j  2  i,j-1  3  i,j 

p.  T*'*  Ax  Ay 

M  .  +  M  v?n}  +  M  vjn? 

=  1  1-1 ,J  2  itJ-1  3  i,J 

(n+1) 

p.  .  Ax  Ay 

1,0 

M  E?n]  .  +  Ml!11!  +  Mjj11! 

=  i  1-1,3.  2  i,J-i _ 

I  m_Ll  1 


More  generally,  if  the  four  sides  of  cell  (i,j)  are  numbered  by 


k  =  1,  2,  3,  and  1  as  shown  below. 


i,j  3 


we  can  define  a  function  C.  .(k)  for  this  cell  referring  to  side  k  such 

i,0 


C.  .(k)  =  1  if  fluid  flows  into  cell  (i,j)  across  side  k 
0 

=  0  if  fluid  flows  out  of  cell  (i,j)  across  side  k 
Then  the  final  velocity  and  energy  are  given  by  the  expressions 

(c.  .0)  u[n]  .  AM<nJ  .  ♦  C.  .(2)  s!n>  ,  AM  <">  , 

1/0  \  1/0  1-1,0  1-2,0  10  1,0-1  1,0-2 


+  c.  .(3)  U.  ,  .  AM,  i  . 

lj  1+1,0  1+2,0 


+  c.  ,(i>)  3<n>  ,  A^n).  ,  +  SU).  -jP!n). 


i,0  i,j+1  i,0+2  i; 


n  n  AxAy  -  [1  -  C  (1)]  (10) 

1  j  J  L  1  j  J  1 


X  AM?n?  .  -  [1  -  C  .(2)]  AM?n?  i 

1-2,0  i,0  1,0-2 


,(3)]  -  [1  -  C.  ,(4)  ]  AmJ.^.JVI 


-  n  -  *i9p»  ^;i,o  - [1  ■  ci,o(4)]  AMi,o+i;;/K,o 'AxAy 


/.nt1}  =  (c.  .(1)  v(.n] 

1,0  y  1,0  1-1, 


.  am!11]  .  + 
0  1-2,0 


c.  .(2)  v!n?  Ah^n^  1 
1,0  1,0-1  1,0-2 


+  c.  .(3)  v(n]  .  Aiin^  + 

1,0  1+1,0  i+lj 


(Equation  continued) 
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C.  .(4)  T<n>  ,  ,  ♦  rtn).  {p'n)  AxAy  - 

1,0  1,0+1  l/J+2  1/0  l  1/0 


[1  -  C,  ,(i)3 


1/0 


X  AM^l  .  -  [1  -  c  (2)3  AM<n)  1  (11) 
:l~2>3  1^J“2 


- 1>  -  4?L  - [i  -  ami”U? 

!i”’: 1  -  (ci,j(,)  AMi-lo  +  Ci,J(2)  AMili 


AxAy 


+  C.  .(3)  E.  .  .A  M^.nl  . 
i,J  i+!/0  1+2/0 


C.  .(4)  E.  .  Am!“!  i  +  E?n).  V.n)  AxAy  -  [1  -  C  (1)3  (1 

i,0  i/J+1  1,0+2  I/O  l  1/0  1/0 


2) 


X  AR^  .  -  [1  -  C.  (2)]  AM^n).  ^ 


-  [1  -  c  4(3)3  AH^}  -  [1  -  c  (4)]AM^n)  1 

ljJ  x+2*o  ±y^  x9 0^2 


( n+1 )  . 

LCi,j  AxAy 


This  completes  the  transition  of  the  fluid  configuration  from  time 
n  to  time  ( n+1 ) . 


C»  Comments 

In  this  section  a  discussion  will  he  given  of  such  questions  as  the 
conservation  of  mass,  energy,  and  momentum  within  the  framework  of  the 
difference  equations,  and  the  stability  and  accuracy  of  the  method. 

1 .  Compliance  with  Conservation  Laws 

Consider  a  rectangular  array  of  cells  covering  a  region  of  fluid 
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under  investigation  and  labeled  by  cell  indices  extending  from  i  =  1  to 

i  =  i  and  j  =  1  to  j  =  j  .  We  want  to  examine  the  change  of  the 
max  max 

total  mass,  energy,  and  momentum  within  the  mesh  between  time  nAt  and 
( n+1 )  A  t .  In  particular,  we  want  to  show  that  the  differencing  scheme 
described  above  is  consistent  with  the  conservation  of  mass,  energy,  and 
momentum. 

a.  Mass  Conservation 

The  change  in  the  mass  of  cell  (i, j)  in  the  transition  from  time 
nAt  to  (n+1 )  At  is 


AM?n}  ■ 


>+1) 


-  P 


(n) 

i>  jJ 


AxAy 


=  Atl"} 


-  A 


+  AM 

If  J“2 


A^n). 


The  mass  change  for  the  entire  system  in 


this  time  interval  is  therefore 


^ - 1 

’LL 


-4-L 


AM^i  . 

1  +  2  JO 


+  AM  ^n)  i 
1fJ~2 


AM 


(n) 

i,  j+i_ 


Since  the  mass  flow  into  cell  (i, j)  from  the  left  is  equal  to  the  flow  out 
of  cell  (i-l,j)  from  the  right,  etc.,  all  terms  in  this  sum  pertaining  to 
boundaries  interior  to  the  mesh  cancel  one  another  leaving  the  result 
^max  -  "jmaxr-  -i 


AM(.n)  x  . 
1max+^,  J 


i=1 


AM^l  -  AM^ 

1J2 


1,  J  +fj 
y  max  *■ 


The  mass  change  of  the  system  in  the  course  of  one  calculation  cycle  is 
seem  to  be  due  to  fluid  entering  or  leaving  across  the  limiting  boundaries 
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of  the  mesh.  No  "creation"  or  "destruction"  of  mass  occurs  in  the  inter¬ 


ior  of  the  mesh. 

b.  Energy  and  Momentum  Conservation 

The  redistribution  of  energy  and  momentum  among  the  mesh  cells  occurs 
in  two  separate  ways  within  a  calculation  cycle:  first  in  Phase  I,  where 
each  cell  is  treated  as  a  Lagrangian  mass  element  and  fluid  is  restrained 
from  flowing,  and  again  in  Phase  III,  where  the  effects  of  the  mass  trans¬ 
port  are  taken  into  account. 

In  Phase  I,  energy  and  momentum  axe  transmitted  to  a  cell  by  means  of 
the  pressure  forces  acting  upon  it.  The  change  in  these  quantities  within 
the  system  is  the  difference  between  the  sum  of  the  tilde  and  the  initial 


energy  and  momenta  of  the  mesh  cells. 


AxAy 


La  La 

i=1  j=1 


AxAy 


Inserting  the  expressions  for  v^.  and  E^  [Egs.  (5,  6,  7)]  into  these 


equations  gives 
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A> 


PT,*'-X 

i,0 


>)  - 
PJ+i  i  pi_i  i 

Ax  Ay  - i-24Li  At 

Ax 


A?T,y'  'I  AxAy^ 


(n)  .  <n) 


As  in  the  case  of  the  mass  flow,  all  terms  pertaining  to  cell  boundaries 
interior  to  the  mesh  appear  twice  —  once  with  a  positive  and  once  with  a 
negative  sign,  and  so  cancel,  tfe  are  left  with  energy  and  momentum  flux 
terms  between  the  boundary  of  the  mesh  and  the  exterior  region  only 
Jmax- 

~  V1  I  (n)  (n)  (n)  (n)  ,  . 

A  E  =  >  Pi  '  ui  '  -  pi  i  ui  i  AyAt 

/L>  L  2>j  2>J  1+2>J  1+2,JJ 


\  (n)  (n)  (n)  (n)  .  -  + 

+  L  Ki  \i  '  pi,j  4vi,J  4  AiAt 

*  *  ’  *  *  max  ■'  max 


.  ~  V1  (n)  (n) 

A  P  =  )  Pi 

x  Aj  L  2,0  i  Q 


AyAt 


A  P  =  V  P(nj  -  p‘n),  AxAt 

y  £j  Li,i  1'Jmax+yJ 

Similarly,  in  Phase  III  mass  crossing  a  cell  boundary  carries  with  it 
an  amount  of  energy  and  momentum  determined  by  the  specific  energy  and  ve¬ 
locity  of  the  cell  from  which  it  is  removed.  As  was  mentioned  in  the  dis¬ 


cussion  of  this  phase  of  the  calculation  cycle,  conservation  of  energy 


22- 


and  momentum  are  used  explicitly  to  ensure  that  energy  and  momentum  lost 
by  one  cell  are  gained  by  another.  (The  equations  that  demonstrate  this 
are  identical  to  those  illustrating  mass  conservation  except  that  now  the 
mass  flows  would  be  multiplied  by  appropriate  specific  energy  or  velocity 
components.)  The  energy  and  momentum  changes  of  the  system  occurring  in 


Phase  III  are 


A  E 


III 


max 


V  UHs<n>  +1^1 

L  L  i,L  1^max+2 

i=1 

4 . 

j=1 


T9*  aM^eC11)  a  1 

AMi  .  E  ,  -  AM.  .J,  a  *v>  ^ 

2^  ^-max^0 


i 

a  p111  =  r  [a  ^  iai  -  ^  ^  td 

x  L  1,2  1,L  1,Jmax+2  1,UJ 


A  P 


III 


max  r  n 


i=1 


a  M(.n|  v<n>  -  aM<n) 

1,2  i>h  1 


n)  ^n)l 

'Vx+^  ^Uj 


(n)  $(“)  _  A  J11)  ,  ^ 

i,3  R/JJ 


+  )  A  M\  v'  '  -  AM'  '  i 

H,  L 

o=1 


where  is  the  tilde  specific  energy  of  the  cell  in  column  i  on 

i,L 
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the  lover  boundary  of  the  mesh  from  vhich  the  mass  is  flawing,  and  simi¬ 
larly  for  the  other  energy  and  velocity  terms.  If  mass  is  flowing  into 
the  system  from  outside  the  mesh,  the  energies  and  velocities  to  be  used 
are  specified  by  the  external  conditions  imposed  on  the  calculation. 

The  total  energy  and  momentum  changes  per  cycle  are  the  sums  of  the 
contributions  to  these  quantities  from  Phases  I  and  III.  Energy  and  mo¬ 
mentum  are  rigorously  conserved  within  the  differencing  scheme  and  their 
changes  can  occur  only  through  interaction  with  the  region  exterior  to 
the  mesh. 

2.  Stability  and  Accuracy 

No  thorough  analysis  of  the  questions  of  the  stability  and  accuracy 
of  the  difference  scheme  described  above  has  been  made.  The  following 
comments  will  therefore  either  be  of  a  self-evident  nature  or  will  be 
based  largely  on  observations  on  past  calculations!  experiments. 

Stability  and  accuracy  are  closely  related.  Both  can  be  affected  by 
relatively  small  changes  in  the  higher  order  derivative  terms  implied  by 
the  replacement  of  the  differential  hydrodynamic  equations  by  difference 
equations.  Basically,  the  difference  scheme  under  consideration  is  accur¬ 
ate  to  first  order  in  the  time  and  space  intervals.  However  true  accur¬ 
acy  at  the  end  of  an  extended  calculation  depends  not  only  on  the  accur¬ 
acy  over  a  single  time  step,  but  on  the  manner  in  which  successive  errors 
tend  to  accumulate  or  cancel.  How  a  specific  form  of  differencing  affects 
this  latter  aspect  of  accuracy  is  not  subject  to  simple  analysis.  In  the 
course  of  developing  the  difference  equations  of  the  previous  section,  a 
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number  of  variations  in  the  form  of  such  quantities  as  the  boundary  pres¬ 
sures  and  the  mass  flows  was  investigated.  The  final  choice  for  these 
quantities  was  based  on  comparisons  of  calculated  and  analytic  solutions 
for  simple  problems  like  the  propagation  and  reflection  of  a  plane  shock 
(see  Fig.  2).  In  this  way,  differencing  procedures  with  large  cumulative 
error  were  hopefully  eliminated. 

A  more  direct  way  in  which  accuracy  is  limited  by  the  magnitudes  of 
the  space  and  time  intervals  is  that  transient  phenomena  occurring  in 
times  shorter  than  At  or  details  of  a  flow  pattern  smaller  than  the  size 
of  a  mesh  cell  will  always  be  washed  out. 

The  difference  equations  are  subject  to  the  Courant  stability  cri- 
k  * 

terion  cAt<  Ax  which  requires  that  a  disturbance  moving  at  the  local 
sound  speed  c  may  only  affect  adjacent  mesh  points  in  a  single  time  cycle. 
A  generally  more  restrictive  condition  for  supersonic  problems  is  that 
u  At  <  Ax  and  vAt  <  Ax  since  the  difference  equations  do  not  provide  a 
mechanism  for  the  transfer  of  fluid  between  nonadjacent  cells. 

Explicit  Eulerian  difference  equations  are  often  unconditionally  un¬ 
stable  with  respect  to  small  perturbations  in  the  flow  quantities,  par¬ 
ticularly  at  low  velocity.  Stability  can  be  achieved  by  the  addition  of 

5 

a  fictitious  viscous  pressure  to  the  hydrostatic  pressure  in  the  momen¬ 
tum  and  energy  equations.  Such  viscous  pressure  terms  cause  a  diffusion 
of  momentum  and  energy  which  has  a  limiting  effect  on  the  gradients  of 
these  quantities  and  constitutes  a  source  of  entropy  production  through 

the  conversion  of  kinetic  to  internal  energy.  Though  the  diffusion 
_ 

In  the  following  we  will  set  Ax  =  Ay. 
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caused  by  these  terms  contributes  to  the  inaccuracy  of  a  difference 
scheme,  it  serves  the  important  function,  in  addition  to  its  stabilizing 
effects,  of  smearing  out  discontinuities  such  as  shock  fronts,  which 

5 

enables  the  calculation  of  these  phenomena  without  special  treatment. 
Typically,  discontinuities  such  as  shock  fronts  are  smeared  over  a 
three-cell  width,  so  that  a  reduction  in  cell  size  will  improve  the  sharp¬ 
ness  of  the  discontinuity. 

The  dissipation  terms  of  the  difference  equations  described  above 
are  implicit  rather  than  explicit  and  arise  from  the  "repartitioning"  of 
the  energy  and  momentum  in  cells  on  the  basis  of  total  energy  and  momen¬ 
tum  conservation  after  the  mass  flow,  i.e.  in  Phase  III  of  the  calcula¬ 
tion.  The  dissipative  mechanisms  involved  can  be  illustrated  by  consider¬ 
ing  a  box  containing  N  distinct  mass  phases  having  different  velocities 
and  internal  energies,  and  then  allowing  these  to  fuse  into  a  single  mass 
with  one  velocity  and  one  internal  energy.  Originally  the  contents  of 
the  box  could  be  described  by  the  set  of  masses,  velocities  and  specific 
internal  energies, 

VVV  k  -  ’’  •••’  H 

The  initial  momentum  and  internal  energy  are 
N 

v* 

k=1 

N 

- 1  “A 

k=1 
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Since  mass,  energy,  and  momentum  must  be  conserved,  the  final  internal 
energy  will  be 


where 


m. 


The  change  in  internal  energy  will  be 


j>k=1 

— *  2  — ►  2  — ►  — »  .  — ♦  —fr  .  2 

Since  vk  +  v^  -  2v^»vk  =  (vk  -  v  J  the  change  in  internal  energy, 

A  IE,  will  be  greater  than  or  equal  to  zero  and  will,  in  fact,  only 
vanish  if  all  N  initial  velocities  are  equal.  Therefore,  the  mixing  of 
the  component  masses  in  general  results  in  a  conversion  of  kinetic  to 
internal  energy,  a  dissipative  entropy  producing  process. 

The  form  of  the  viscous  pressure  implicit  in  the  difference  equa¬ 
tions  may  be  obtained  by  converting  these  back  to  differential  equa¬ 
tions  at  the  point  (i, j),1  retaining  terms  through  first  order  in  Ax 
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and  Ay.  That  is,  we  let 


pi,0  -,p’ 


em  "E’  ui,j  -u>  vi,j 


and  replace  flow  quantities  at  neighboring  points  by  appropriate  Taylor 
expansions  about  (i, j) .  In  performing  this  reduction  it  is  convenient 
to  choose  a  definite  flow  direction,  say  in  the  positive  u  and  v  direc¬ 
tions.  The  resultant  differential  equations  are 


+  div(pu}  =  0 

+  +  It  =  I;  (l  puAxtO  +  b  % pvAy^) 

+  div(pv3  +  ^  =  | -  (^  Pv^y^)  +  p^x  §*) 

+  div(pEu!)  +  div  (pi?)  =  (l  puAx^)  +  ^  PvAy^) 

Apart  from  the  first  order  terms  on  the  right,  these  are  the  Eulerian 
hydrodynamic  equations.  That  they  can  be  obtained  from  the  difference 
equations  justifies  the  three  phase  difference  scheme. 

The  viscous  pressure  is  seen  from  the  momentum  equations  to  have 
the  appearance  of  a  tensorial  pressure  resulting  from  viscous  interactions, 

* 

The  first  order  terms  in  At  do  not  seem  to  have  as  important  a  role  as 
the  viscosity  terms,  undoubtedly  because  At  «  Ax.  They  are  also  pri¬ 
marily  diffusion- like  quantities. 
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1 

a  du 

A  ^ 

TAy^ 

2P 

Sv 

bv 

vAy^  _ 

However,  written  in  the  diadic  form  —  pu  Ar*Vu,  where  Ar  =  Axi 

•4 

+  Ayj  ,  it  is  seen  that  this  pressure  is  not  an  invariant  form  due  to  the 
presence  of  the  fixed  vector,  Ar!  Its  elements  are  large  only  when  the 
velocity  and  its  gradients  are  large.  (The  diffusion  terms  in  the 
energy  equation  above  Eire  also  small  when  v  is  small. )  Thus  in  regions 
of  near  stagnation  the  viscosity  will  be  ineffective  and  the  instabili¬ 
ties  of  the  difference  equations  will  cause  an  exponential  growth  of  any 
perturbations  to  the  solution.  However,  once  the  velocity  has  reached  a 
magnitude  comparable  to  the  local  sound  speed,  the  viscosity  will  act  to 
limit  further  growth  of  the  instability.  This  type  of  bounded  instabili¬ 
ty  is  illustrated  by  Fig.  2a  where  the  calculated  velocity  for  a  one  di¬ 
mensional  steady  shock  reflected  from  a  barrier  is  shown.  The  shock 
front  is  moving  to  the  left  and  the  velocity  behind  it  should  be  zero. 
However,  instabilities  cause  large  oscillations  in  this  quantity. 

The  difficulties  with  the  difference  equations  at  low  velocities  can 
be  alleviated  by  supplementing  the  implicit  viscous  pressure  with  an  ex¬ 
plicit  term  in  the  momentum  and  energy  equations.  A  satisfactory  way  of 
accomplishing  this  is  to  add  such  a  term  to  the  pressure  at  a  cell  bound¬ 
ary  if  the  fluid  speed  is  found  to  be  small,  compared  to  the  sound  speed 
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at  this  boundary.  For  example,  assuming  the  fluid  is  a  perfect  gas  and 
considering  the  right-hand  boundary  of  cell  (i, j),  one  could,  make  a  test 
to  see  if 


20Li ,j  =  *7  -  ’>  (li,3  +  +  Vi,/  +  (Ti,j  +  Vi,/1 


03) 


where  c  1  .is  the  sound  speed  at  the  boundary  and  A  is  a  constant  grea- 
i+2>«J 

ter  than  one.  If  the  inequality  is  satisfied  no  additional  viscous  term 

6 

is  necessary.  If  it  is  not,  a  Landshoff  type  viscous  pressure  could  be 
added  to  the  boundary  pressure,  i.e.. 


(pi,J  +  +  pi+1,^"Bci+|,j  **  +  Vl.j1  ‘  V1,J* 

w 

With  such  a  prescription  the  difference  equations  will  remain  stable  at 
fluid  speeds  of  a  few  percent  of  the  local  speeds  of  sound.  Optimum  con¬ 
stants  A  and  B  can  be  determined  from  tests  with  simple  one  dimensional 
problems.  For  the  calculations  in  Section  IV  of  this  report  where  the 
fluid  was  air  and  y  -  1  .4,  A  and  B  were  1 .5  and  0.3,  respectively. 

Figures  2b  and  2a  compare  the  calculated  velocity  of  a  steady  shock  which 
has  been  reflected  from  a  rigid  wall,  with  and  without  a  supplementary 
viscous  pressure.  The  Improved  stability  is  self  evident. 

As  a  final  remark,  it  might  be  noted  that  negative  internal  energies 
can  arise  under  special  circumstances  in  the  course  of  a  calculation  as 
in  the  PIC  method.  This  is  particularly  true  during  the  initial  time 
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X 


Figure  2a.  Plot  of  a  calculated  velocity  profile  vs  position  for  a  plane  shock  reflected 
from  a  rigid  vail  using  no  explicit  viscous  pressure.  The  true  shock  position 
is  indicated  by  a  vertical  line  vith  an  arrow  specifying  direction  of  motion. 


A1I0013A 


■52. 


as  in  Fig.  2a,  but  with  an  explicit  viscous  pressure 


cycles  when  a  cold  fluid  at  rest  is  acted  upon  by  some  driving  force. 


In  this  case,  the  source  of  the  negative  internal  energy  can  be  found 

in  Eqs.  (5),  (6),  and  (7)  of  the  Phase  I  calculation.  If  initially 

u(n)  v(n)  and  are  zero  but  a  pressure  gradient  exists,  u^  and 

i>j  i>J 


v^  will  not  vanish  but  will  in  general  remain  zero  since  the 

i,0 

time  n  velocities  are  used  to  calculate  the  energy  fluxes.  Therefore 
one  will  obtain 


since  =  0. 

ij  j 


This  result  is  due  to  a  somewhat  improper  phasing  he- 


tween  the  momentum  Eqs.  (5)  and  (6)  and  the  energy  Eq.  (7)  and  could 


be  largely  corrected  by  replacing  the  velocities  in  the  latter  by  the 
just  calculated  tilde  velocities.  However  these  negative  internal  ener¬ 
gies  do  not  seem  to  persist  for  more  than  a  few  time  cycles  and  have  not 
led  to  any  difficulties,  so  their  occasional  occurrence  has  been  ignored. 


III.  BOUNDARY  CONDITIONS  AND  THE  TREATMENT  OP  SEVERAL  FLUIDS 

For  a  single  fluid  the  most  common  types  of  boundary  conditions 
fall  into  two  categories,  those  at  open  ends  of  the  mesh  across  which 
fluid  can  enter  or  leave  and  those  at  rigid  barriers  within  or  at  the 
extremities  of  the  mesh.  The  different  cases  which  fall  under  each  cate¬ 
gory  will  be  enumerated  below. 
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A.  Boundary  Conditions  at  Open  Ends  of  a  Mesh 
1 .  Prescribed  Input 

Suppose  fluid  is  flowing  into  a  rectangular  mesh  from  the  left-hand 
side.  Then  the  left  side  (^,j)  of  cell  (1,  j)  forms  the  boundary  with 
the  exterior  region.  It  is  seen  from  Eqs.  (5),  (7)*  (8),  (9)#  and  ( 10)  that 
the  following  quantities  must  be  known  in  order  to  apply  the  difference 
equations: 


Pi  *t  v!  A  , 

2fJ  2 to 


That  is,  we  must  know  the  energy  and  momentum  flux  into  the  cell  from  the 
outside  and  the  mass  transport  quantities.  The  indices  (0,  j)  and  (-1 , j) 
refer  to  fictitious  cells  to  the  left  of  mesh  cell  (1,j).  The  input 
boundary  conditions  are  prescribed  by  specifying  the  above  eight  quanti¬ 
ties  for  each  time  cycle.  More  simply,  it  is  usually  possible  to  pre¬ 
scribe  the  nature  of  the  flow  to  the  left  of  cell  ( 1  ,j) .  If  p_-j  y 
v  .,  p  ,  u  ,  v  ,  E  are  given  for  the  fictitious  external  cells 

“  *  j  J 

(-1,j)  and  (0,j)  each  cycle,  the  boundary  quantities  at  (^,j)  can  be  cal¬ 
culated  with  the  equations  of  Section  II. 


2.  Continuous  Outflow 

In  order  to  allow  fluid  to  flow  out  of  a  mesh,  the  external  boundary 
must  be  sufficiently  far  from  any  source  of  disturbance  so  that  the  flow 
pattern  can  be  extrapolated  to  the  outside.  The  simplest  example  of  this 
situation  is  where  the  flow  is  uniform,  say  at  the  right-hand  end  of  the 
mesh.  If  the  mesh  is  N  cells  long,  the  conditions  in  a  fictitious  cell 


(N  +  1 ,  j)  to  the  right  of  the  last  cell  (N,  j)  should  always  be  the  same 
as  in  (N,j).  Suitable  boundary  conditions  can  be  achieved  by  allowing 
for  cell  (N  +  l,j)  and  setting 


\,y  VN+1,j  VN,o' 


E  =  E 

n+1,j  N,j 


at  the  end  of  Phase  I  of  each  cycle  of  calculation,  and  setting 


(n+1)  =  (n+1) 
pN+1,j  pN,j 


J>+1)  =  (n+1) 


(n+1) 

N+1, 


.(n+1) 

’N+1, 


=  E 


(n+1) 

N,j 


at  the  end  of  Phase  III. 

B.  Interactions  with  Rigid  Barriers 
1 .  Barriers  at  Cell  Boundaries 

Barriers  at  cell  boundaries  occur  most  commonly  when  the  mesh  con¬ 
sists  of  a  channel  to  which  the  fluid  is  confined  or  in  cylindrical  co¬ 
ordinates  where  the  axis  of  symmetry  can  be  represented  by  a  rigid  bar¬ 
rier.  The  boundary  conditions  are  determined  by  the  necessity  that  there 
must  be  no  energy  flux  or  fluid  flow  across  the  barrier.  This  will  be 
assured  if  the  velocity  normal  to  the  barrier  is  zero.  The  momentum  flux 

at  the  obstacle  boundary  is  not  zero  and  can  be  obtained  approximately 
by  extrapolating  the  pressure  to  the  wall.  If  the  mesh  is  fine,  the 
pressure  at  the  barrier  can  be  chosen  as  that  at  the  center  of  the  cell 
bounded  by  the  obstacle. 

By  way  of  an  example,  assume  the  right-hand  side  of  cell  (i,j)  re¬ 
presents  the  boundary  with  a  rigid  obstacle.  As  in  the  case  of  the 
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open-ended  boundary,  a  fictitious  mesh  cell  (i  +  1 ,  j)  can  be  assumed  and. 
its  state  prescribed  in  terms  of  that  of  cell  (i,j)  so  that  application 
of  the  equation  of  Section  II  will  automatically  produce  the  correct 
boundary  conditions.  The  prescription  for  cell  (i  +  1 ,  j)  is  that  the 
density  and  specific  energy  must  at  all  times  be  the  same  as  in  cell 
(i,j),  but  the  velocity  must  be  the  mirror  image  of  that  in  (i, j) .  Thus 
at  the  end  of  Phase  I  one  can  set 


V  =  v 

i+l,J 


and  at  the  end  of  Phase  III, 


P 


(n+1) 

1+1,  J 


-  E 

1+1,  J 


(n+1) 

i,J 


9 


U 


(n+1)  = 

1+1,  J 


u(n+1) 

i,J  1 


v(n+1)  =  (n+1) 
1+1,3  i,j 


In  the  case  that  the  obstacle  has  a  corner  as  shown  by  the  shaded  area 


below. 

i+l,j+l 

'.j 

777777Z 

W///A 

W/ 

the  fictitious  cell  (i  +  1,3)  will  be  assigned  two  different  sets  of 
variables  depending  on  whether  computations  are  being  done  in  cell  (1,3) 
or  ( i  +  1 ,  j+1). 

The  use  of  a  fictitious  cell  with  assigned  variables  is  not  actually 

necessary  though  it  often  is  convenient.  In  practice  one  can  test  to 

see  if  calculations  are  being  made  in  a  "boundary  cell"  like  (  i,3)  shown 

~(n)  ~(n+1 )  _  _ 

above,  and  if  so,  merely  set  p.^j  =  and  u.^  =  u.+^.  -  0. 
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If  this  second  procedure  is  followed,  it  has  been  found  that  setting 

u  =  +  u  for  use  in  the  mass  flow  calculation  improves  the  ac- 

i+1,j  i,J 

curacy  of  the  computations. 

2.  Arbitrary  Linear  Barriers  Within  a  Cell 

As  in  the  previous  case  of  a  barrier  at  a  cell  boundary,  the  dif¬ 
ference  equations  applicable  when  the  barrier  is  within  a  cell  are  de¬ 
termined  by  the  requirements  that  no  energy  flux  or  mass  may  cross  the 
boundary.  In  fact,  treatment  of  a  barrier  at  a  cell  boundary  is  merely 

a  special  case  of  what  will  be  discussed  here. 

* 

Placing  a  barrier  within  a  mesh  cell  has  two  primary  effects.  It 
shifts  the  center  of  mass  from  the  cell  center  to  some  point  nearer  the 
cell  boundary,  and  it  reduces  the  effective  length  of  the  cell  in  at 
least  one  direction.  In  general,  the  hydrodynamic  parameters  describing 
the  state  of  the  fluid  within  a  mesh  cell  refer  specifically  to  the 
center  of  mass  of  the  cell.  When  making  an  interpolation  of  some  quan¬ 
tity  between  cells,  one  in  fact  should  interpolate  between  centers  of 
mass.  For  non- fractional  cells  the  center  of  mass  will  either  coincide 
with  the  cell  center,  as  in  Cartesian  coordinates,  or  be  sufficiently 
close  to  the  cell  center  so  that  the  difference  can  be  neglected.  An 
example  of  the  latter  situation  is  cylindrical  coordinates  where  even  for 
cells  adjacent  to  the  axis  of  symmetry  the  center  of  mass  differs  from 
the  cell  center  by  only  two  tenths  of  the  radial  cell  dimension.  Since 

Cells  containing  or  adjacent  to  a  barrier  will  be  called  "fractional 
cells." 
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the  center  of  mass  of  fractional  cells  are  in  general  veil  displaced 
from  the  cell  center,  we  should  by  rights  calculate  those  quantities 
needed  at  cell  boundaries  by  interpolation  with  respect  to  centers  of 
mass.  However  the  errors  involved  in  interpolating  with  respect  to  cell 
centers  are  felt  to  be  minor,  even  for  fractional  cells;  and  so  in  prac¬ 
tice,  these  effects  of  center-of-mass  displacement  have  been  neglected. 

The  reduction  of  the  effective  cell  dimensions  in  fractional  cells 
is  a  more  serious  matter,  particularly  if  the  center-of-mass  displace¬ 
ments  axe  disregarded.  This  is  because  the  stability  requirements  for 
the  difference  equations  at  such  a  cell  may  now  be  violated  with  the  re¬ 
sultant  appearance  of  large  localized  fluctuations.  Therefore  special, 
procedures  must  be  adopted  to  minimize  these  difficulties. 

There  are  twelve  distinct  types  of  fractional  cells  containing  sin¬ 
gle  linear  barriers  which  must  be  treated  individually.  They  are  of  the 
following  form  (Fig.  3),  where  the  shaded  area  represents  the  presence 
of  a  rigid  obstacle: 


(I)  (2)  (3)  (4) 


15)  (6)  (7)  (8) 


(9)  (10)  (II)  (12) 


Figure  3.  Types  of  fractional  cells. 


For  each  fractional  cell  a  set  of  five  geometrical  quantities  f. 

J 

A.  !  ,  A.  i,  A.  !  ,  and  A  .  i  will  be  needed,  f  .  is  the  fraction 

of  the  total  cell  "volume,"  Ax  Ay,  which  is  occupied  by  fluid;  A^i^  is 
the  fraction  of  the  "area"  of  cell  side  (i  -  ir,j)  which  is  open  to  flow; 


etc. 


Using  the  above  quantities,  difference  equations  for  fractional 
cells  can  be  obtained  from  those  for  "full"  cells  given  in  Section  II  B 
by  multiplying  mass  and  energy  transport  terms  at  a  given  cell  side  by 
the  corresponding  fractional  area.  A,  and  by  multiplying  the  densities 
occurring  explicitly  in  Eqs.  (5*  6,  7>  8,  10,  11,  12)  by  f.  .•  The  pres- 
sures  in  Eqs.  (5)  and  (6)  now  represent  averages  of  pressures  at  the  perim¬ 
eter  of  the  region  occupied  by  fluid  in  the  fractional  cell.  To  avoid 


overly  complicating  the  computing  procedure,  pressures  have  been  calcu¬ 
lated  as  if  the  fractional  cell  were  effectively  rectangular,  with  the 
pressure  at  cell  boundaries  taken,  as  before,  to  be  the  simple  average  of 
the  pressures  in  the  cells  on  either  side,  and  the  pressure  at  the  obsta¬ 
cle  face  being  taken  as  that  within  the  fractional  cell.  The  pressure 
difference  in  Eq.  (5)  must  be  multiplied  by  the  larger  of  the  fractional 
areas  A.  1  ,  and  A.  x  ,,  and  that  of  Eq.  (6)  by  the  larger  of  A  A  and 

Ai>j+2* 

With  these  changes,  Eqs.  (5)  through  (9)  become 

=  ^(n)  _  Pi+?.j  ~  MaX(Ai+^J,  (5t) 


Ax 


7=) 
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(n)  „(n) 
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AMi+l  1 
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In  Eqs.  (10),  (11),  and  (12)  we  merely  replace  each  p.  .  by  f  .  p.  ..  All 
other  equations,  for  instance  those  defining  the  Pjj£|.  j  and  ui+f,  j  aPP6^” 
ing  in  Eq.  (9)  are  treated  in  the  same  manner  as  in  the  case  of  full  cells 
with  a  barrier  along  a  cell  boundary,  as  described  in  III  B1  above. 

Another  procedure  which  can  be  used  without  undue  error  in  the  calcu¬ 
lation  of  tilde  quantities  in  fractional  cells  is  to  regard  the  cell  as 

a/  a/ 

full-fractional  as  in  III  Bl  for  the  purpose  of  obtaining  u^  and  v±^. 
This  produces  a  somewhat  different  transfer  of  momentum  to  the  obstacle 
than  results  from  the  use  of  Eqs.  (5’)  and  (6')#  but  does  not  seriously 
affect  the  flow  pattern  if  the  mesh  is  not  too  coarse.  This  alternate 
differencing  cannot  be  applied  to  the  calculation  of  E.  .  in  fractional 
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cells  since  it  would  result  in  a  violation  of  energy  conservation.  No 
energy  can  be  transferred  to  the  obstacle. 

In  some  instances  the  fractional  volumes  of  cells  of  type  (1),  (2), 
(5),  or  (4)  may  be  sufficiently  small  so  that  stability  requirements  are 
violated.  The  second  treatment  of  the  tilde  velocities  is  then  particu¬ 
larly  useful.  In  addition,  that  for  the  tilde  energy  must  be  changed 
from  Eq.  (7')  to  prevent  instabilities.  One  such  modification,  which  is 
not  in  contradiction  with  energy  conservation,  is  to  allow  only  the  frac¬ 
tion  f  .  of  the  energy  crossing  the  fully  open  boundary  of  these  types 
J 

of  cells  to  enter  or  leave  the  fractional  cell,  the  remaining  fraction 
(1  -  f  .)  being  returned  to  the  cell  from  which  the  energy  flowed.  For 
example,  for  a  type  (5)  cell  shown  below  the  E  equations  would  become: 


The  calculation  of  a  viscous  pressure  at  the  boundaries  of  fractional 
cells  can  be  accomplished  similarly  to  that  for  a  normal  cell  with  the 
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following  additional  prescription: 


a.  If  the  cell  side  has  part  of  its  area  open  for  fluid  flow  into 
an  adjacent  cell,  calculate  the  viscous  pressure  terms  as  if 
this  side  were  entirely  unobstructed. 

b.  If  a  side  of  cell  (i,j)  is  entirely  closed  to  fluid  flow,  Eqs. 
(13)  and  (14)  should  be  replaced  by 


12  ^  A 

2  °i,j  <  A 


2 

+  Vi,J- 


and 

(pi,j  +pi+i,j) 


Ti,j 


(u1) 


(IV) 


The  treatment  of  fractional  cells  described  above  is  rather  crude  and 
can  be  improved  upon  in  many  ways.  However  it  has  the  virtue  of  relative 
simplicity  and  does  not  deter  from  accuracy  as  long  as  the  Eulerian  mesh 
is  made  fine  compared  to  the  obstacle. 


C.  Two  Material  Cells 

The  presence  of  several  distinct  fluids  separated  by  interfaces  within 
the  confines  of  an  Eulerian  mesh  increases  the  difficulty  of  such  calcula¬ 
tions  many  fold.  Special  treatment,  similar  to  that  for  "fractional"  cells 
in  one  fluid  problems,  must  be  given  to  cells  containing  more  than  one 
fluid  at  each  phase  of  the  calculation;  and,  in  addition,  a  method  for  mov¬ 
ing  the  interface  must  be  prescribed.  In  the  following,  a  method  for  per¬ 
forming  Eulerian  calculations  in  cells  containing  two  fluids  will  be  out¬ 
lined.  This  procedure  has  been  tested  with  reasonable  success  in  several 
trial  problems.  It  is  characterized  by  an  indirect  method  for  moving  the 
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fluid  interface.  This  is  as  opposed  to  a  direct  movement  of  a  boundary 
by  means  of  the  explicit  motion  of  points  which  specify  its  position. 

A  single  fluid  in  an  Euler i an  mesh  cell  is  characterized  by  a  den¬ 
sity,  velocity,  and  internal  energy  "located"  at  the  cell  center.  When 
two  fluids  are  present  in  such  a  cell,  each  will  have  a  distinct  density, 
internal  energy,  and  parameter  specifying  the  fraction  of  the  cell  each 
occupies.  A  distinct  velocity  may  also  be  given  to  each  material,  with 
the  restriction  that  the  components  normal  to  the  interface  be  equal, 
in  order  to  allow  for  slippage  at  the  interface.  This  situation  will 
not  be  considered  here,  and  a  single  velocity  for  both  fluids  in  a  given 
cell  will  be  assumed.  In  addition  the  assumption  of  continuity  of  pres¬ 
sure  will  be  made  so  that  there  is  a  single  pressure  in  the  cell  possessed 
by  both  fluids.  The  velocity  and  pressure  of  a  multi-fluid  cell  is  thus 
treated  as  in  a  single  fluid  cell  and  is  assumed  "located"  at  the  cell 
center.  Although  this  point  will  in  general  no  longer  be  the  cell  center 
of  mass,  errors  due  to  this  deviation  should  not  be  important «  Finally, 
the  interface  segment  within  a  single  cell  will  always  be  assumed  to  be 
a  straight  line.  A  specification  of  which  pair  of  cell  boundaries  this 
line  intersects  and  the  position  of  the  intersection  points  are  required 
to  complete  the  information  about  the  two  fluid  cell. 

The  complexity  of  the  equations  of  motion  for  a  two  fluid  cell  depend 
on  the  degree  of  precision  one  wants  to  incorporate  in  them.  Ideally  each 
material  in  a  two  fluid  cell  should  be  treated  independently  as  a  single 
fluid  in  a  "fractional  cell"  except  for  such  constraints  as  continuity  of 


pressure.  This  will  in  fact  he  done  in  Phases  II  and  III  of  the  calcu¬ 
lation,  hut  the  cell  will  he  treated,  here  as  a  single  full  cell  in 
Phase  I. 

Consider  an  arbitrary  two  fluid  cell  such  as  that  shown  below: 

— J-Ai+I/2,J(,) 

_ ^  A 1+1/2,  )(2) 

As  with  a  single  fluid  cell,  it  is  assumed  that  at  time  nAt  the  den¬ 
sity,  mass,  and  specific  energy  of  both  fluids  are  known  and  that  a 
single  pressure  and  velocity  has  also  been  assigned  to  each  cell.  The 
Phase  I  equations,  in  which  no  mass  flow  occurs,  can  be  deduced  from 


Eqs.  (2)  and  (3)  to  be 
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where  f^n^(1)  and  are  the  fractional  volumes  of  cell  (i, j)  at 

1/0  1/0 

"time"  n  occupied  by  fluids  (1)  and  (2),  respectively.  If  R^n^0 )  and 

R(n)(2)  are  the  respective  masses  of  the  two  fluids  in  cell  (i, j)  divided 
1/0 


by  the  volume  AxAy,  then 
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and 


The  R  . ,  and  p  . ,  are  not  all  independent  since 
i fO  s  i.»J  s 

f$n).0)  +  f(,n),(2)  =  1 

Equation  (17)  does  not  permit  an  independent  determination  of  E, 

E^n^(2) .  Unless  these  are  calculated  in  a  manner  similar  to  that  for 
fractional  one  fluid  cells  discussed  above,  an  additional  assumption  must 
be  made  at  this  point.  One  such  is  that  the  energy  brought  into  the  cell 
by  the  flux  terms  in  Eq.  (17)  is  distributed  uniformly  across  the  cell. 
Then  each  fluid  is  incremented  in  energy  by  an  amount  proportional  to  its 
fractional  volume  so  that 
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a  =  1,2 


The  assumption  leading  to  the  above  equation  is,  of  course,  arbitrary. 

The  calculation  of  the  mass  flows  in  Phase  II  is  similar  to  that  for 
the  single  fluid  fractional  cell  described  above.  There  are  now  two  sep¬ 
arate  sets  of  mass  flows,  however,  one  for  each  fluid.  In  order  to  avoid 
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extrapolation  of  densities  beyond  fluid  interfaces,  mass  flows  in  multi- 
fluid  problems  have  used  boundary  densities  which  are  simply  the  density 
of  the  cell  from  which  the  fluid  flows.  Talcing  into  account  the  area 
available  for  each  fluid  flow  across  a  given  cell  boundary  as  determined  by 
the  position  of  the  interface  in  the  cell,  we  have  for  instance  for  side 


(i+l,j)  of  the  two  fluid  cell  pictured  above: 
If 


and 
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a  =  1,2 


and  otherwise 


In  these  equations,  A^Ui  .(a)  is  the  fractional  area  of  side  (i+2>  j) 

i+zfj 

of  cell  (i,j)  available  for  the  flow  of  fluid  a.  For  two  fluids 

4:L<’>  +A<i%Js)  =  1 


If  a  fluid  has  no  contact  with  a  cell  boundary,  its  fractional  area  at 
that  boundary  will  be  zero. 

These  considerations  of  the  mass  flow  across  side  (i+§,j)  of  the 
cell  illustrated  above  can  easily  be  extended  to  the  other  three  sides 
and  to  other  two  fluid  configurations.  The  procedure  can  be  seen  to  have 
one  serious  limitation.  If  there  were  flow  upward  across  the  top, 

(i,j+g),  °f  cell  (i,j),  all  of  mass  (l)  would  have  to  vacate  this  cell 
before  any  of  mass  (2)  could  enter  cell  (i,j+1) .  Thus  the  method  of  cal¬ 
culation  introduces  local  irregularities  into  -the  interface  surface,  though 
on  the  average  the  boundary  will  follow  its  proper  course.  This  difficulty 
can  be  remedied  by  calculating  the  maos  flows  on  the  basis  of  the  volume  of 
one  cell  which  is  transported  to  another.  The  mass  flow  would  then  just 
be  equal  to  the  masses  originally  contained  in  this  small  volume. 

The  final  step  in  the  mass  flow  calculation  for  a  two  fluid  cell  is 
to  check  that  the  existing  masses  do  not  exceed  the  total  mass  original  ly 
in  the  cell.  If  it  should  for  one  of  the  fluids,  each  outgoing  mass  flow 
for  that  material  must  be  reduced  proportionally  to  its  originally  calcu¬ 
lated  value  so  that  the  sum  of  the  outgoing  flows  exactly  equals  the  total 
mass  originally  present.  For  example,  if  the  total  mass  of  material  (1)  is 
M.  .0)  and  the  outgoing  mass  flows  are 
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then  the  replacement  must  he  made 


AMi,j-i0)  -  BA*W,) 


AM.,  .(1)  -4  rAM  i  ,(1) 

i+2,j  i+2 ti 


where 


r  a 


V’1 


AMi4,Ju'-A,W) 


The  third  phase  of  the  two  fluid  calculation  is  essentially  the  same 
as  in  the  single  material  case.  Equations  (8),  (10),  (11),  and  (12)  are 
still  relevant  except  that  it  is  the  "mass"  (per  volume  Ax  Ay) 

n(n)  _  -(n)  (n)  (19\ 

Ri,j  -  fi,j  pi,j 

which  must  replace  the  densities  p.  ,  which  explicitly  appear  in  these 
equations.  It  should  he  noted  that  for  cells  containing  only  a  single  ma¬ 
terial  \  while  for  multi -material  cells  the  new  densities 

p(n+1  gxe  still  unknown  at  this  point  of  the  calculation.  The  third 

phase  of  the  calculation  must  further  contain  a  test,  after  the  calculation 

of  the  new  R  quantities,  to  see  whether  a  given  cell  has  remained  a  one 
J 

of  two  fluid  cell,  or  whether  a  single  material  cell  has  become  a  two 
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material  cell,  or  vice  versa.  If  no  change  in  cell  structure  has  oc¬ 
curred,  no  further  work  is  necessary.  But  if  a  one  material  cell  has 
become  a  two  material  cell,  provision  must  be  made  to  determine  the 
"type"  of  cell  in  accordance  with  Fig.  3,  locate  the  interface,  and  cal¬ 
culate  the  densities;  while  if  the  reverse  is  true,  these  special  con¬ 
siderations  are  no  longer  necessary. 

As  was  stated  above,  at  the  end  of  the  Phase  III  calculation  the 
densities  in  multifluid  cells  are  still  unknown  and  must  be  obtained. 

In  addition  the  new  position  of  the  boundary  between  materials  must  be 
found  and  a  method  given  for  computing  the  pressure  in  these  cells. 

These  quantities  could  be  calculated  directly  on  the  basis  of  a  procedure 
which  involves  the  explicit  movement  of  fluid  interface  with  interpolated 
velocities  during  the  Phase  II  calculations.  However  the  necessity  of 
keeping  the  interface  motion  compatible  with  the  mass  flows  causes  this 
to  be  a  complicated,  though  not  intrinsically  difficult  process.  As  an 
alternative,  one  can  first  obtain  the  pressure  and  densities  in  the  cell 
as  a  consequence  of  the  assumption  of  continuity  of  pressure  and  then  use 
this  information  to  locate  the  intercepts  of  the  interfluid  boundary  with 
the  cell  sides.  By  way  of  exanple,  if  the  two  fluids  were  perfect  gases 
obeying  the  equations  of  state, 

P(1)  =  (^  -  1)  pO)  1(1) 

and 

P(2)  =  b2  -  1)  P(2)  1(2) 

where  I  is  the  specific  internal  energy,  the  continuity  of  pressure  in  the 
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two  material  cell  (i>tj)  requires 


Here  all  quantities  but  the  p.  .(a)  are  known  at  the  end  of  the  Phase  III 

1,  0 

calculation.  Using  the  relations 
R,  Act) 


and 


pi,j(a)  °  q  j(a)  ' 


fi,J0)  +  fi,J(2)  ‘  ’ 


a  *  1,2 


one  can  obtain  an  equation  for  f.  .0) 

R,  ,0) 


R,  .(2) 


<7,  -  ’>  fHrr  v0  ‘ (ra  • 11  i--^7n 


whose  solution  is 


f.  -(1)  =  7 - 

i,0  (7, 


(r,  -  D  Rii;10)  ilit)0> 

-  1)  R<  JT7  t4  fn  +  lr'  -  1)  R,  Jzl  irp) 


1,0 


i,0 


i,0 


Once  f  (1)  is  known,  p.  .0),  p.  .(2),  and  p.  can  be  found  immediately. 
i,j  i,0  i,0  i,0 

In  practice,  if  the  equations  of  state  for  the  materials  are  compli¬ 
cated  and  the  cell  size  small  compared  to  the  distance  over  which  signifi¬ 
cant  changes  in  parameters  occur,  it  may  be  adequate  and  easier  to  extra¬ 
polate  the  densities  in  multifluid  cells  from  those  in  single  fluid  cells. 
The  cell  pressure  can  then  be  obtained  from  those  of  the  two  materials  as 

an  average  weighted  by  the  fractional  volumes,  f.  (ct) ,  of  each. 

-  i,0 

With  the  fractional  volumes,  f.  .(a),  of  the  two  fluids  in  the 

i,0 
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interface  cells  known,  a  simple  procedure  for  obtaining  the  interface 
position  can  be  prescribed  provided  the  slope  of  this  boundary  does  not 
vary  rapidly  from  cell  to  cell.  If  this  latter  condition  is  not  met, 
the  mesh  is  probably  too  coarse  to  provide  an  accurate  calculation  in 
any  event.  This  procedure,  which  will  be  described  below  by  an  example, 
can  be  used  in  Cartesian  coordinates  or  in  cylindrical  coordinates  away 
from  the  symmetry  axis.  Consider  the  situation  illustrated  below: 


FLUID  2 

(M,j) 

bAx 

o+i,j> 

^aAy{^ 

FLUID  1 

We  want  to  calculate  the  intercept  distances  yAy  on  the  right  of  cell 
(i,j).  It  is  assumed  that  the  cell  types,  i.e.,  the  sides  cut  by  the 
interface  in  each  cell,  have  already  been  determined  by  testing  the  com¬ 
position  of  neighboring  cells.  A  "distance"  y  can  be  obtained  in  three 
different  ways  by  the  assumptions: 

(1)  the  interface  in  cell  (i-1,j-1)  and  ( i- 1 , j )  is  a  single  straight  line 

(2)  the  interface  in  cell  (i-1,j)  and  (i,j)  is  a  single  straight  line 

(3)  the  interface  in  cell  (i,j)  and  (i+1,j)  is  a  single  straight  line 


Under  the  first  assumption  one  obtains  the  equations, 

(1)  =?yi  (1-b) 
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a  yT 


b  1  -  b 


Therefore 


*  ■ 2  [w,J  ♦VW,) 

Similarly  under  assumptions  (2)  and  (3)  one  gets 


-  2fe7TT 


y3  =  21 


3  f,  4(1)  -  f,  .1  «0) 


A  simple  average  of  these  three  quantities  will  generally  give  an  ade¬ 
quate  approximate  value  for  the  intercept  distance  y.  Under  some  cir¬ 
cumstances  the  fractional  volumes  f .  .  may  not  be  compatible  with  reason- 

J 

able  straight  line  interfaces  through  a  pair  of  cells.  If  this  is  the 
case,  the  final  value  for  the  boundary  position  may  be  negative  or  grea¬ 
ter  than  1 .  If  this  should  happen,  a  cruder  approximation  for  y  above 
can  be  used  which  is  not  subject  to  this  type  of  error.  This  approxima¬ 
tion  is  to  make  the  boundary  position  equal  to  half  the  sum  of  the  frac¬ 
tional  volumes  of  the  two  adjacent  cells  whose  common  side  is  cut  by  the 
interface.  For  instance  in  the  above  example,  set 


-iff 

2  Li-1, 


0)  +  f.  <0) 


After  the  computation  of  the  intersections  of  the  cell  boundaries 
with  the  fluid  interface,  which  gives  the  areas  (in  Cartesian  coordi¬ 
nates)  available  for  flow  in  bi-fluid  cells,  the  cycle  of  calculation  is 
completed . 
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Appendix  I  contains  a  table  of  equations  for  interface  positions  for 
different  material  configurations  under  the  assumption  that  the  boundary 
though  two  adjacent  cells  is  a  single  straight  line. 

IV.  CALCULATIONS  OF  SUPERSONIC  FLOW  PAST  BLUNT  CYLINDRICAL  OBSTACLES 
As  an  illustration  of  the  numerical  procedures  described  above  and  a 
check  on  their  accuracy,  a  code  for  the  IBM  7090  digital  computer  was 
written  to  calculate  the  flow  past  arbitrarily  shaped  cylindrical  objects. 
This  is  a  problem  which  is  not  generally  amenable  to  direct  analytic  calcu¬ 
lations  and  is  of  some  intrinsic  interest. 

Because  of  the  assumed  axial  symmetry  of  the  obstacles,  cylindrical 
coordinates  have  been  used  for  this  calculation  requiring  several  modifi¬ 
cations  of  the  one  fluid  Eulerian  difference  equations  discussed  in  Sec¬ 
tion  II.  These  changes  are  due  to  the  increase  of  volume  and  effective 
boundary  area  of  a  cell  as  its  distance  from  the  axis  increases.  Equa¬ 
tions  referring  strictly  to  motion  in  the  axial  direction  ( z-direction) 
are  identical  to  the  Cartesian  equations. 

If  i  and  j  are  the  cell  indices  referring  to  the  z  and  r  directions, 
respectively,  and  if  j  =  1  gives  the  row  closest  to  the  symmetry  axis,  the 
complete  set  of  difference  equations  for  non- fractional  cells  in  cylindri¬ 
cal  coordinates  is  as  follows: 

Phase  I 


radial  velocity 


longitudinal  velocity 

(n)  (n) 

?(n)  »  v(n)  _  Pi+p.J  ~  At 


i,j  1,0 


A  z 


r^r 

Ji  >0 


specific  energy 


E 


<■) ,  .  \>  -iV -1:U :  -S-* 

1,0  i,j  L  U  -  i)Ar 


Jn)  v(n)  _  v(n! 

1  l+?,j  1“2,  J _ iliEisL 

A  z 


At 


°i,0 


where 


u(nl  .  =  i  fp?n).  +  V{n]  .1,  etc. 

- i+2> 0  2  L  i;J  1+1,0 J 


Phase  II 


A  M 


(n) 

i+i,  0 


A  »4n>  x 
1,0  2 


(j 
(  J 
(j 


-  4!L  tL  - r2  - 


- 1>  4“L  i%i  a  r2  At 

-  1)  P$n),  i  i  ArAzAt 


a  4nLi =  j  4nL  4ni+i ArA  zAt 

1,0+2  l,<3+2  1,0+2 


where  the  p^l  .,  v^l  etc.  are  calculated  as  in  Section  II. 
1"2>0  1“2,J 
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Phase  III 


(n+1)  =  (n)  1 _ 

*'J  (j4)Ar2Az 


L4r-St 


♦  m'n),  a  -  A  M<“)  1  . 

I,  j-2 


4"1’  -  (v”  %i,3  +  ci,j(2)  AMS-* 


+  C  (5)  X  AmH 


+  Ci,j(4)^4  +  *[*]  {(J4)  Ar2Azp[“j  -  P-C^O)]  AM^i 


-  [1-C  (2)]AMj.n)  i  -  [1-C  (3)]  AMJJ 

1#0  J”2  1>J  1+2>J 


Xn) 


[,-ci,/4)1  AM$J/ 


/  .  i\.  2  (n+1 ) 

(0-5)Ar  Azp) 

J-tJ 


where  X  is  u,  v,  or  Ej  and  C.  ,(k)  is  one  or  zero  depending  on 

j 

whether  fluid  is  flowing  into  or  out  of  cell  (i,j)  across  side 


k  as  defined  in  Section  II. 


The  differences  between  the  Cartesian  and  cylindrical  equations  are 
seen  to  be  minor.  Explicit  artificial  viscosities  can  be  calculated  as 
explained  earlier  and  the  rules  for  treating  fractional  cells  apply  un¬ 
changed  except  for  the  form  of  the  fractional  volumes  and  areas.  Appen¬ 
dix  II  contains  a  table  of  fractional  areas  and  volumes  for  the  different 
types  of  fractional  cells  in  cylindrical  coordinates.  One  additional  dif¬ 
ference  that  should  be  noted  is  the  boundary  condition  at  the  axis  of 
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symmetry.  This  must  be  treated  as  a  rigid  barrier. 

The  problems  which  were  run  with  the  obstacle  code  were  chosen  to 

2 

coincide  with  a  set  of  supersonic  wind  tunnel  experiments  by  Marschner 

3 

and  a  similar  experiment  at  a  higher  Mach  number  by  Oguchi.  The 
Marschner  experiments  which  were  simulated  were  the  flow  of  air  past  a 
right  circular  cylinder  at  Mach  1 .58  and  the  flow  of  air  past  a  cone  of 
75°  vertex  angle  attached  to  a  circular  cylinder,  at  Mach  numbers  1 .41 , 
I.58,  and  1.71.  The  Mach  1.41  case  was  also  rerun  at  double  the  mesh 
resolution  in  order  to  obtain  an  idea  of  gain  in  accuracy  which  could  be 
brought  about  by  decreasing  the  cell  size.  The  experiment  of  Oguchi 
which  was  considered  was  again  the  flow  of  air  past  a  right  circular  cyl¬ 
inder  but  at  Mach  4.3.  In  al 1  of  these  experiments,  a  detached  shock 
wave  was  present.  The  data  on  the  flows  which  were  available  and  with 
which  the  machine  computations  were  compared  were  the  positions  of  the 
shock  fronts  and  the  pressure  distributions  at  the  obstacle.  Oguchi' s 
results  also  provide  the  position  of  the  sonic  line  —  line  of  Mach  1  — 
and  the  stream  line  pattern  immediately  beyond  the  shock. 

The  hydrodynamic  differencing  scheme  which  was  described  in  the 
previous  sections  is  explicitly  a  method  for  time  dependent  problems. 

For  stationary  problems  it  is  necessary  to  begin  the  calculation  with  an 
arbitrary  configuration  of  the  fluid,  always  keeping  the  correct  boundary 
conditions,  and  advance  the  time  until  a  steady  state  is  reached.  Thus 
the  computation  contains  the  dynamics  of  the  approach  to  steady  state  as 
well  as  the  steady-state  flow  itself.  In  the  present  calculations  a  rec¬ 
tangular  mesh  was  set  up  as  shown  with  the  axis  of  symmetry  as  the 
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bottom  boundary.  The  upper  boundary  was  assumed  rigid. 


Fluid  entered  at  the  left  with  the  free  stream  conditions 

o=1,  u  =  0,  c  =  sound  speed  =  1 , 
p0  *  0  ’  0  * 

v  =  M„  =  Mach  number 
0  0 

and  excited  at  the  right.  The  equation  of  state  was  that  for  a  polytropi 
gas  with  y  =  1.4. 

p  =  (7  -  1)  p[E  -  (u2  +  v2)] 

The  initial  fluid  parameters  throughout  the  mesh  were  chosen  equal  to  the 
free  stream  parameters.  These  initial  conditions  correspond  to  the  impul 
sive  start  up  of  the  obstacle  in  a  stationary  medium  as  seen  from  a  re¬ 
ference  frame  in  which  the  obstacle  is  at  rest.  In  the  case  of  the  cone 
tipped  cylinders,  velocities  throughout  the  mesh  were  sufficiently  high 
that  the  use  of  an  explicit  viscous  pressure  was  not  considered  essential 
Such  a  term  was  used,  however,  with  the  blunt  nosed(l80°  cone)  cylinders. 
The  following  table  lists  the  mesh  characteristics  and  viscous  pressure 
parameters  used  in  the  present  calculations: 
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Mesh  Size  Obstacle  Viscous  Pressure 


Cone  Angle  Mach  Number  (r  x  z) 

Radius 

Parameters 

75° 

1.41 

40  X  60  cells 

9  cells 

none 

75° 

1.58 

40  x  60  cells 

9  cells 

none 

75° 

1.71 

40  x  60  cells 

9  cells 

none 

75° 

180° 

1.41 

1 .58 

60  X  60  cells 

44  x  44  cells 

18  cells 

14  cells 

none 

A  =  1 .5 

B  =  0.5 

180° 

4.5 

44  x  44  cells 

14  cells 

A  =  1.5 

B  =  0.5 

The  time 

steps  At  were 

chosen  so  that  M  X  At 

~  1/10  Ar, 

In  all  cases 

Ar  and  Az  were  equal. 

Problems  were  allowed  to  run  until 

it  was  judged 

that  changes  in  the  mesh  parameters  between  a  time  t  and  t  +  40  X  At 
were  negligible.  These  running  times  were  on  the  order  of  two  to  three 
hours. 

In  order  to  compare  the  results  of  the  machine  calculations  with  the 
experimental  data  for  each  problem,  plots  were  made  of  the  computed 
position  of  the  detached  shocks  and  of  the  ratio  of  the  pressure  at  the 
surface  of  the  obstacle  to  the  true  stagnation  pressure.  In  addition, 
stream  lines  for  each  flow  were  obtained  by  following  the  trajectories 
of  "test  particles"  through  the  mesh  using  the  computed  steady-state 
velocities.  Linear  interpolation  between  the  velocities  at  the  cell 
centers  was  used  to  obtain  the  test  particle  velocity.  The  stream- line 
plots  which  follow  were  made  directly  from  computer  results  by  an  auto¬ 
matic  plotter  and  not  hand  drawn.  In  the  Mach  4.5  case,  an  additional 
graph  compares  the  computed  and  experimental  stream  lines. 
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In  the  case  of  the  blunt  cylinders,  surface  pressures  were  obtained 
by  a  linear  extrapolation  from  the  pressures  at  the  centers  of  the  cells 
immediately  in  front  of  the  obstacle.  A  similar  procedure  was  followed  for 
the  75°  cones.  However  there,  in  an  attempt  to  minimize  the  effects  of 
fluctuations  of  pressure  in  the  small  fractional  cells  and  to  provide  a  sys¬ 
tematic  method  for  extrapolation  to  the  obstacle,  two  rows  of  cells  were 
constructed  parallel  to  the  surface  of  the  cone  as  shown  below.  Pressures 


at  the  centers  of  the  tilted  cells  were  obtained  from  the  pressures  within 
the  normal  mesh  as  an  average  weighted  by  the  area  of  each  normal  cell 
within  a  slant  cell.  With  these  averaged  values,  the  pressures  along  the 
cone  were  obtained  by  extrapolation.  In  the  graphs,  distance  along  the 
cone  is  measured  in  -units  of  the  slant  length  s  as  shown.  This  length  is 
just  the  radius  of  the  cylinder  for  the  blunt  faced  obstacles.  The  solid 
lines  give  the  experimental  curves  while  the  calculated  values  are  shown 
by  circled  points  connected  by  broken  lines.  As  mentioned  before,  pres¬ 
sures  were  plotted  in  units  of  the  true  stagnation  pressure  in  each  case. 
These  stagnation  pressures  can  be  obtained  by  following  an  element  of 
fluid  along  the  axial  stream  line.  From  the  extreme  left  until  the  shock 
front,  this  element  flows  with  the  free  stream  parameters.  At  the  shock 
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the  parameters  change  discontinuously  according  to  the  normal  Hugoniot 
relations  to 


2  +  (7-I )Mq 


p  V 

i-00 

1  9 


v1  = 


p’  =  po  +  povo(vo  •  T,) 


From  just  beyond  the  shock  to  the  stagnation  point,  v  =  0,  Bernoulli ' s 
law  and  the  adiabatic  equation  of  state  may  be  used  to  give  the  stagna¬ 
tion  pressure  and  density.  Thus 


7 


7-1 


p' 

p' 


and 


or 


The  stagnation  pressures  for  the  situations  considered  are  as  follows: 


Mach  Number 

Stagnation  Pressure 

1.41 

2.20 

1 .58 

2.66 

1.71 

5.05 

4.5 

17.34 

Shock  fronts  were  obtained  from  the  calculations  by  first  making 


plots  of  density  versus  axial  position  for  each  row  of  cells  in  the  case 

o 

of  the  blunt  cylinders  and  Mach  number  versus  axial  position  for  the  75 
cones.  These  show  a  steep  rise  in  the  vicinity  of  the  shock  followed  by 
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a  more  gradual  rise  or  decline.  The  plots  are  given  in  Appendix  III.  Be¬ 
cause  of  the  smearing  of  the  shock  front  by  the  implicit  viscous  pressure 
inherent  in  the  calculation,  the  position  of  the  shock  is  necessarily  sub¬ 
jective  to  some  degree.  The  shock  front  lies  roughly  halfway  between  the 
point  of  the  initial  steep  rise  and  the  point  where  the  rise  is  judged  to 
terminate  in  the  density  or  Mach  number  curve.  A  set  of  unconnected  points 
marking  the  predicted  shock  position  have  been  plotted  in  each  case  as  well 
as  a  solid  line  marking  the  experimental  shock  front.  In  addition,  a  pair 
of  dashed  lines  obtained  from  the  approximate  one  quarter  and  three  quar¬ 
ter  positions  on  the  calculated  shock  rise  curves  have  been  shown.  These 
define  an  approximate  limiting  region  in  which  the  predicted  shock  must 
lie.  In  the  Mach  4.3  problem,  the  experimental  and  theoretical  sonic  lines 
are  also  given.  Figures  4  through  22  show  the  stream  lines,  shock  front, 
and  surface  pressure  distribution  for  each  of  the  six  problems  considered. 

An  examination  of  the  graphs  shows  that,  in  general,  reasonably  good 
agreement  between  theory  and  experiment  can  be  obtained  by  the  numerical 
methods  employed.  The  estimated  positions  of  the  shock  fronts  agree  with 
experiment  to  better  than  one  cell  width  in  all  cases  except  that  of  the 
blunt  cylinder  at  Mach  1.58.  There  the  predicted  shock  lies  farther  from 
the  obstacle  than  that  determined  experimentally,  though  both  curves  are 
similar.  The  displacement  of  the  calculated  shock  intercept  with  the 
axis  of  symmetry  with  respect  to  the  experimental,  value  amounts  to  about 

(Text  continues  Page  81 ) 
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EXPERIMENTAL  STREAM  LINES 


COMPUTED  STREAM  LINES 


ftp^ 

/  /  /  / 
/  /  / 

/ 


-SHOCK  FRONT 


Figure  20.  Stream  lines,  180  cone.  Mach  4.3 


z/R 


Figure  22.  Detached  shock  position,  I80°eone.  Mach  4.3 
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a  five  percent  deviation.  The  source  of  this  discrepancy  is  not  known, 
especially  since  no  estimates  of  experimental  errors  were  given  by 
Marschner.2 

The  comparison  of  the  calculated  and  experimental  pressure  distribu¬ 
tions  at  the  obstacle  is  somewhat  marred  by  the  relatively  large  fluctua¬ 
tions  for  the  75  cones  of  nine  cell  radius.  Increasing  the  resolution 
of  the  mesh  as  was  done  in  the  case  of  the  Mach  1 .41  cone  reduced  these 
fluctuations  substantially  and  resulted  in  good  agreement  with  experi¬ 
ment.  In  general,  the  calculated  pressures  follow  the  trend  of  the  ex¬ 
perimental  pressures  at  roughly  the  correct  magnitudes.  The  use  of  an 
explicit  viscous  pressure  term  in  the  cone  problems  would  have  probably 
reduced  the  pressure  fluctuations  near  the  obstacle  and  generally  im¬ 
proved  the  agreement  with  experiment.  For  the  blunt  cylinders,  no  ex¬ 
perimental  pressures  at  Mach  1 .58  were  available.  In  the  Mach  4.5  prob¬ 
lem  agreement  between  calculation  and  experiment  is  good  except  near  the 

shoulder  of  the  cylinder  where  the  calculated  pressures  are  somewhat 
high. 

In  conclusion,  the  present  numerical  method  for  solving  the  hydro- 
dynamic  equations  for  supersonic  compressible  flows  will  give  adequate 
results  where  extreme  precision  is  not  required.  The  accuracy  of  the 
scheme  depends  on  the  nature  of  the  problem  and  on  specific  quantity 
under  consideration  as  well  as  on  the  relative  size  of  the  mesh  compared 
to  the  expected  spatial  variations  of  physical  quantities,  on  the  time 
step,  and  on  the  proper  use  of  an  explicit  viscous  pressure  term,  in 
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the  present  problems,  this  accuracy  may  be  estimated  at  about  five  per¬ 
cent.  There  is  still  much  room  for  improvement  in  the  method.  Higher 
order  accuracy  could  probably  be  incorporated  in  the  equations  without 
undue  difficulty.  However,  perhaps  the  most  important  undertaking  would 
be  to  improve  the  treatment  of  boundary  condition,  particularly  in  the 
case  of  fractional  and  multifluid  cells* 
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APPENDIX  I 


The  following  equations  give  the  intercept  with  the  cell  boundaries 
of  the  interface  between  a  pair  of  fluids  in  an  Eulerian  mesh  under  the 

following  assumptions: 

a.  Cartesian  coordinates  are  appropriate. 

b.  The  cell  dimensions  are  Ax  =  Ay  =  1 . 

c.  The  fluid  interface  is  a  single  straight  line  segment  within  a 
pair  of  adjacent  mesh  cells. 

d.  The  fractional  areas  f .  .  occupied  by  each  fluid  within  a  two 

i,0 

fluid  cell  are  known.  The  subscripts  (i,j)  will  be  abbreviated 
by  the  single  capital  letters  A  or  B. 


fA0)  =  \  (a  +  b) 
fg(  1  )  “  5  0>  +  c) 

b  =  £  (a  +  c) 
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2  1f*0)  +  fJ1) 


r-'-i 


APPENDIX  II 

The  following  equations  give  the  fractional  areas,  A^,  and  frac¬ 
tional  volumes,  f,  available  to  a  fluid  in  a  mesh  cell  partially  occupied 
by  a  rigid  obstacle  in  cylindrical  coordinates.  It  is  assumed  that 
Ar=  Az  =  1  and  that  the  left-  and  right-hand  cell  sides  will  be  desig¬ 
nated  by  subscripts  1  and  J,  and  the  bottom  and  top  sides  by  2  and  4, 
respectively. 
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zo0-r  )  [j  -  1/3  (1-r  )] 


11. 


APPENDIX  III 

The  following  six  graphs  show  the  changes  in  the  Mach  number  of  den¬ 
sity  in  the  region  of  the  detached  shock  for  the  situations  considered  in 
this  report.  On  each  graph  a  series  of  plots  are  given  for  the  appro¬ 
priate  quantity  in  different  mesh  rows.  The  position  z  equals  zero  re¬ 
presents  the  left-most  extremity  of  the  obstacle  in  the  mesh.  The  esti¬ 
mated  shock  position  is  marked  by  an  x  on  each  plot  and  the  approximate 
bounding  region  in  which  the  shock  front  should  lie  is  marked  by  small 


circles 
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Figure  III.1 .  75  cone,  Mach  1.41,  40  x  60  mesh 


Figure  III. 4.  75  cone,  Mach  1.41,  60  X  60  mesh 


DENSITY  CURVES  ARE  DRAWN  FOR 


80  cone,  Mach  1 


FHE  DENSITY  CURVES  ARE  DRAWN  FOR 


Figure  III. 6.  180  cone,  Mach  4.5 


